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ABSTRACT 


An  adaptive  atmospheric  flow  model  is  described  and  results  of  integrations  with  this  model  are  presented. 
The  adaptive  technique  employed  is  that  of  Berger  and  Oliger.  The  technique  uses  a  finite  difference  method 
to  integrate  the  dynamical  equations  first  on  a  coarse  grid  and  then  on  finer  grids  which  have  been  placed 
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based  on  a  Richardson-type  estimate  of  the  truncation  error  in  the  coarse  grid  solution.  By  correctly 
coupling  the  integrations  on  the  various  grids,  periodically  re-estimating  the  error  and  recreating  the  finer 
grids,  uniformily  accurate  solutions  are  economically  produced.  The  '"primitive^ hydrostatic  equations  of 
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meteorology  are  solved  for  the  advection  of  a  barotropic  cyclone  and  for  the  development  of  a  baroclinic 


disturbance  which  results  from  the  perturbation  of  an  unstable  jet.  These  integrations  demonstrate  the 


feasibility  of  using  multiple,  rotated,  overlapping  fine  grids.  Direct  computations  of  the  truncation  error  are 
used  to  confirm  the  accuracy  of  the  Richardson-type  truncation  error  estimates. 


This  work  has  been  supported  by  the  Office  of  Naval  Research  under  contracts  N00014-84-K-0267  and 
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1.  INTRODUCTION 

Accurately  computing  atmospheric  flows  is  a  difficult  task.  Integrations  must 
be  performed  accurately  and  in  a  timely  manner  if  the  simulations  are  to  be  helpful 
to  weather  forecasters,  but  the  variety  of  scales  which  must  be  included  for  successful 
representation  of  the  flow,  and  about  which  forecasters  need  information,  strain  the 
resolution  capabilities  and  speed  of  the  most  advanced  computers.  Planetary  waves 
are  many  thousands  of  kilometers  long  and  exist  for  weeks.  At  the  opposite  end  of 
the  scale  of  meterological  interest  are  the  turbulent  eddies  which  have  length  scales 
of  a  few  centimeters  and  last  only  seconds.  We  axe  far  from  resolving  all  scales 
in  a  single  computation.  Typically  we  resolve  only  a  small  subset  of  the  scales  of 
motion,  parameterize  others  and  disregard  what  is  left.  Unfortunately  this  is  often 
not  adequate. 

There  are  several  approaches  that  are  presently  being  used  to  address  the  res¬ 
olution  problem  in  numerical  weather  prediction  (NWP).  Grids  are  often  vertically 
nonuniform,  with  more  layers  close  to  the  surface.  Different  vertical  coordinates 
have  been  used  to  achieve  this  in  a  more  natural  way.  For  example,  the  poten¬ 
tial  temperature  (B  =  T(paurface/p)R/cf )  may  be  used  as  the  vertical  coordinate 
with  a  resulting  increase  in  resolution  around  fronts.  The  normalized  pressure 
o  =  pj Psurface  is  more  commonly  used  as  a  vertical  coordinate.  This  system  al¬ 
lows  increased  accuracy  in  the  specification  of  the  surface  boundary  because  of  the 
simplicity  of  the  surface  boundary  condition.  Horizontal  layers  smoothly  follow  the 
surface.  These  approaches  are  passive  methods  for  increasing  resolution.  They  are 
a  direct  result  of  the  form  of  the  equations  or  difference  methods. 

Another  technique  for  increasing  the  resolution  in  atmospheric  prediction  codes 
consists  of  placing  finer  grids  inside  a  coarser  grid  at  any  location  where  greater 
resolution  or  accuracy  is  desired.  In  the  atmospheric  science  literature  the  fine 
grids  are  known  as  nested  grids  and  implementation  has  been  achieved  in  two  forms, 
one-way  interactive  and  two  way  interactive. 

One  way  interaction  is  the  simplest  nested  grid  approach.  The  lateral  boundary 
conditions  for  the  fine  grid  are  supplied  from  the  coarse  grid  solution.  Information 
is  passed  from  the  coarse  grid  to  the  fine  grid  but  not  from  fine  to  coarse,  hence 
it  is  called  one-way  interactive.  Many  regional  atmospheric  models  are  one-way 
interactive  for  they  receive  their  lateral  boundary  values  from  global  model  solutions 
but  have  no  effect  on  these  solutions.  The  larger  scales  of  the  flow  which  cannot  be 
simulated  on  the  fine  grid  are  allowed  to  affect  the  fine  grid  solution.  The  major 


V.V.VAA 


•V  A'  v-'k  v\,  v.a  ,vv 4v.Vjy 


problem  is  the  different  wave  speeds  which  result  from  the  different  resolutions  on 
the  two  grids.  Discontinuities  and  distortions  can  develop  at  the  fine  grid  boundary 
as  a  result  of  having  no  feedback  from  the  fine  to  the  coarse  grid  (Haltiner  and 
Williams,  1980).  Sponge-type  boundary  conditions  have  been  developed  for  one¬ 
way  interactive  fine  grid  models  which  make  possible  useful  results. 

The  one-way  interactive  approach  implicitly  assumes  that  small  scale  phenom¬ 
ena  have  no  major  influence  on  the  larger  scale  flow  treated  on  the  coarse  grid.  This 
is  not  generally  true  for  there  are  two-way  exchanges  of  energy  between  scales.  The 
two-way  interactive  nested  grid  approach  addresses  this  problem.  The  procedure 
consists  of  integrating  the  fine  grid  along  with  the  coarse  grid.  The  lateral  boundary 
conditions  for  the  fine  grid  are  taken  from  the  coarse  grid  solution.  The  solution 
on  the  coarse  gnd  is  updated  with  the  fine  grid  solution  at  any  location  where  the 
two  coincide.  The  scheme  is  thus  two-way  interactive  because  the  fine  grid  solution 
does  influence  the  coarse  grid  solution. 

The  two-way  interactive  method  has  been  implemented  in  two  forms,  one  where 
the  fine  grid  location  is  fixed  and  another  where  the  fine  grid  is  allowed  to  move 
during  the  time  integration.  An  example  of  the  former  is  the  Nested  Grid  Model 
(NGM)  which  has  been  delevoped  at  the  National  Meteorological  Center.  The  Na¬ 
tional  Weather  Service  distributes  NGM  results  as  guidance  to  forecasters.  The 
model  consists  of  a  hemispherical  grid  over  which  two  finer  grids  have  been  placed. 
The  fine  grids  are  centered  over  North  America  because  this  is  where  forecast  in¬ 
formation  is  needed. 

Several  tropical  cyclone  models  have  fine  grids  which  move  during  the  integra¬ 
tion  to  keep  the  fine  grid  over  the  cyclone.  This  is  accomplished  by  moving  the  fine 
grid  when  a  solution  feature,  such  as  the  surface  low  associated  with  the  cyclone, 
moves.  In  all  cases  the  fine  grids  are  aligned  with  the  coarse  grids,  but  they  may 
move  incrementally  up,  down  or  sideways.  Examples  of  these  are  tropical  cyclone 
models  described  by  Harrison  (1973)  and  Jones  (1977). 

The  tropical  cyclone  models  attempt  to  provide  resolution  around  features 
which  are  local  and  poorly  resolved  by  the  coarse  grid.  The  NGM  provides  increased 
resolution  over  a  continent  because  increased  resolution  may  be  necessary  there, 
but  it  is  difficult  to  know,  a  priori,  precisely  where  it  will  be  necessary.  From 
a  global  perspective  tropical  cyclones,  fronts,  jet  streams  and  other  atmospheric 
phenomena  are  spatially  and  temporally  localized  and  are  often  poorly  resolved  in 
present  atmospheric  models.  Local  phenomena  should  be  handled  adaptively  but 
no  adaptive  atmospheric  flow  solvers  exist.  The  nested  tropical  cyclone  models  are 
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not  truly  adaptive.  The  initial  location  of  the  cyclone  must  be  known  to  locate  the 
fine  grid.  The  fine  grid  will  remain  even  if  the  cyclone  disap  pears,  and  if  a  new 
cyclone  were  to  appear  elsewhere  no  new  fine  grid  would  appear  over  it. 

Phenomena  in  many  other  fluid  flows  which  are  difficult  to  resolve  are  often 
localized.  Adaptive  solvers  do  exist  for  many  of  these  flows.  In  general  two  adaptive 
strategies  axe  used.  In  the  first  all  existing  gridpoints  axe  redistributed  from  regions 
of  small  solution  variation  to  regions  of  large  solution  variation.  These  global  meth¬ 
ods  vary  in  the  criteria  and  methods  used  to  move  the  points,  but  in  all  cases  the 
total  number  of  points  remains  the  same.  They  axe  usually  used  in  conjunction  with 
grid  transformation  methods  which  involve  mapping  an  irregular  physical  domain 
into  a  rectangular  computational  domain.  The  second  strategy  involves  adding  or 
deleting  grid  points  so  as  to  obtain  a  desired  solution  accuracy.  The  additions  and 
deletions  are  local,  thus  the  techniques  are  local  grid  refinement  techniques. 

Atmospheric  flows  appear  ideally  suited  to  local  grid  refinement  techniques  be¬ 
cause  the  important  phenomena  are  localized.  The  local  grid  refinement  techniques 
can  be  broken  into  two  catagories:  One  in  which  the  new  points  are  inserted  or 
imbedded  into  the  existing  grid,  and  hence  only  one  grid  exists,  and  a  second  where 
refinements  are  placed  over  the  existing  grid,  the  refinements  constituting  separate 
grids. 

A  fine  example  of  embedding  new  points  on  an  existing  grid  is  the  work  of 
Dannenhoffer  and  Baron  (1986).  Their  code  solves  transonic  flow  over  a  2-D  air¬ 
foil.  Refinements  are  based  on  refinement  parameters,  for  example  first  or  second 
order  differences  in  the  density,  pressure  or  entropy.  An  expert  system  handles  the 
refinement  parameters  and  rules  governing  how  and  where  to  refine. 

In  this  technique  grids  are  no  longer  rectangular  in  nature  and  neither  is  the 
data  structure  which  holds  the  solution  fields.  A  significant  amount  of  information 
must  be  stored  to  describe  the  grid  structure.  The  solver  which  uses  this  grid 
structure  is  complex.  Even  with  this  complexity  and  loss  of  rectangulaxity  much 
vectorization  of  the  code  is  possible  and  efficient  integrations  are  being  performed. 

An  example  where  refinements  are  placed  over  the  existing  grid  is  the  scheme  of 
Berger  and  Oliger  (1984).  The  same  scheme  has  been  used  by  Berger  and  Jameson 
(1985)  who  also  solve  transonic  flow  over  a  2-D  airfoil.  In  this  technique  the  re¬ 
finements  axe  sepaxate  rectangular  grids  rather  than  being  points  embedded  in  the 
coarse  grid.  Any  solver  which  works  on  a  rectangle  can  be  used,  because  the  solver 
is  just  a  module  called  by  the  adaptive  routines  to  advance  the  solution.  Berger 
and  Jameson  calculate  results  similar  to  those  of  Dannenhoffer  and  Baron.  There 
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are  many  other  differences  between  the  two  schemes  and  interested  readers  should 
consult  the  referenced  papers. 

Solver  complexity  is  a  strong  barrier  in  NWP  to  techniques  which  imbed  points 
into  existing  grids.  Weather  prediction  codes  solve  much  more  than  a  simple  set 
of  dynamical  equations.  There  are  equations  for  water  vapor  (and  possibly  water 
in  its  other  states),  routines  which  calculate  radiational  heating  and  heating  due  to 
phase  changes  of  water,  routines  which  model  cloud  effects,  complex  calculations 
for  fluxes  of  heat  and  moisture  into  the  atmosphere,  and  usually  parameterizations 
of  other  processes.  Most  codes  are  the  result  of  many  peoples’  effort  over  several 
years.  Proven  and  tested  routines  are  often  borrowed  from  one  code  to  put  into 
another.  Adaptive  methods  using  refinements  which  are  separate  grids  appear  the 
logical  choice  for  use  as  the  basis  of  an  adaptive  weather  model.  Existing,  well 
tested  software  can  often  be  used  with  only  minor  changes  and  procedures  can  be 
written  with  little  knowledge  of  the  adaptive  routines. 

In  this  paper  we  present  results  from  an  adaptive  atmospheric  flow  solver  which 
uses  the  method  and  software  developed  by  Berger  and  Oliger  (1984).  The  adap¬ 
tive  method  operates  on  multiple,  component  grids.  Fine  grids,  which  overlie  the 
coarse  grid(s),  are  created  and  removed  based  on  a  Richardson-type  estimate  of  the 
truncation  error  in  the  finite  difference  solution.  The  goal  is  to  maintain  a  given 
accuracy  for  a  minimum  amount  of  work.  This  is  the  first  attempt  with  this  method 
to  adaptively  solve  a  system  describing  a  three-dimensional  time-  dependent  flow. 

Our  purpose  is  to  show  that  an  adaptive  atmospheric  flow  solver  is  feasible  and 
that  the  adaptive  technique  will  produce  self-consistent  results.  In  essence  we  are 
proving  a  concept,  the  concept  being  (1)  that  refinement  should  occur  only  where 
necessary,  as  dictated  by  the  error  in  the  numerical  solution,  and  in  this  way  improve 
the  accuracy  and  overall  resolution  of  the  entire  solution  and  (2)  that  this  can  be 
accomplished  by  using  the  method  of  Berger  and  Oliger.  Hence,  we  wish  only  to 
demonstrate  that  our  adaptive  model  is  successful  when  compared  with  results  from 
the  same  solver  on  a  single  grid,  this  being  sufficient  to  demonstrate  the  feasibility 
of  the  adaptive  atmospheric  flow  solver.  For  prediction  purposes  one  would  use  the 
best  available  solver  for  the  scales  one  is  attempting  to  predict. 

In  section  2  we  review  the  adaptive  grid  refinement  technique  of  Berger  and 
Oliger.  The  solution  procedure  is  outlined  along  with  a  description  of  the  data 
structure  and  program  design.  Section  3  describes  the  equations  and  solver  we  have 
used  in  the  adaptive  package.  The  initial  test  cases  describe  flow  in  an  idealized 
atmosphere:  adiabati'  flow  in  a  periodic  channel  with  no  moisture  present  in  the 
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atmosphere  and  no  diabatic  heating.  We  are  using  the  Euler  equations  with  the 
hydrostatic  assumption,  the  so  called  “primitive  equations”  of  meteorology.  At  the 
end  of  the  section  we  discuss  the  issues  of  stability  and  accuracy  for  these  equations 
as  they  are  used  in  this  adaptive  context.  In  section  4  we  examine  the  results  of  two 
simulations,  one  for  a  barotropic  cyclone  and  another  for  a  baroclinically  unstable 
jet,  and  show  that  our  adaptive  model  is  self  consistent  and  successful  in  simulating 
these  flows.  We  examine  the  error  estimate  in  section  5  and  conclusions  follow  in 
section  6. 


2.  REVIEW  OF  THE  ADAPTIVE  GRID  REFINEMENT 
TECHNIQUE 

We  describe  the  adaptive  procedure  as  used  for  2-D  hyperbolic  problems.  For 
large  scale  atmospheric  flows  the  horizontal  scales  axe  orders  of  magnitude  larger 
than  the  vertical  scales.  In  our  adaptive  model  we  refine  in  the  horizontal  and  keep 
the  same  number  of  layers  in  the  vertical.  Thus  we  treat  atmospheric  flow  as  a  2-D 
grid  refinement  problem  even  though  it  is  a  3-D  flow.  We  are  currently  working  on 
a  model  which  will  also  refine  in  the  vertical. 

The  adaptive  method  is  based  on  the  idea  of  using  multiple,  component  grids 
on  which  the  partial  differential  equations  axe  solved.  Refined  grids  Eire  created 
or  removed  based  on  a  Richardson-type  estimate  of  the  truncation  error  in  the 
finite  difference  solution.  The  goad  is  to  maintain  a  given  accuracy  in  the  numerical 
solution  for  a  minimum  amount  of  work.  A  complete  description  of  the  method  can 
be  found  in  Berger  (1982)  and  Berger  and  Oliger  (1984). 

The  solution  procedure  for  the  adaptive  grid  method  is  as  follows.  We  begin 
with  a  solution  on  a  coarse  grid  that  has  been  integrated  to  some  time  t.  The  error 
introduced  through  the  use  of  the  numerical  procedure  is  estimated  at  grid  points 
and  where  these  errors  are  judged  to  be  too  large  the  points  are  flagged.  Then  2-D 
rectangular  grids  with  finer  meshes  are  fitted  around  these  flagged  points.  These 
subgrids  may  have  orientations  differing  from  the  coarse  grid.  Initial  and  boundary 
conditions  for  these  new  fine  grids  are  taken  from  the  coarse  grid  solution  and  the 
fine  grids  are  integrated  along  with  the  coarse  grid  to  a  new  time  t  +  Atc.  where 
A<c  is  the  coarse  grid  time  step.  Smaller  time  steps  are  taken  on  the  fine  grids 
to  keep  Ax/ At  constant.  The  coarse  grid  solution  is  then  updated  with  the  more 
accurate  fine  grid  solutions.  Errors  may  also  be  estimated  on  the  fine  grids  and  still 
finer  grids  introduced.  Thus,  there  can  be  several  levels  of  fine  grids.  The  errors  on 
the  grids  can  be  re-estimated  every  few  time  steps  and  then  new  fine  grids  can  be 
created  and  old  fine  grids  removed. 

The  large  errors  in  the  numerical  solution  are  usually  associated  with  sharp 
variations  in  the  solutions,  e.g.,  at  fronts  or  other  disturbances.  By  re-estimating 
the  error  after  several  time  steps  and  regridding  we  create  a  scheme  whereby  the  fine 
grids  track  the  disturbances.  The  fine  grids'  arbitrary  orientation  allows  them  to 
line  up  with  the  disturbances  which  minimizes  the  size  of  the  fine  grids  and  provides 
for  better  resolution  in  numerical  schemes. 

The  program  can  be  viewed  as  consisting  of  three  components;  1)  a  data  struc- 
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ture,  2)  a  solver  and  3)  management  routines.  Due  to  the  constructs  of  FORTRAN 
the  data  structure  is  fixed.  It  stores  the  solution  vectors  for  all  grids  and  informa¬ 
tion  about  these  grids.  The  information  and  where  it  is  stored  are  altered  by  the 
management  routines.  These  routines  also  pass  to  the  solver  the  locations  of  the 
solution  vectors  in  the  data  structure.  The  overall  program  can  be  viewed  as  the 
interaction  between  the  solver  and  data  structure  with  the  management  routines 
controlling  this  interaction  and  performing  the  necessary  intergrid  communication 
(setting  boundary  conditions  and  updating). 

The  key  algorithms  axe  those  which  perform  the  integration,  error  estimation 
and  grid  generation.  To  illustrate  how  they  work,  consider  the  grid  arrangement 
shown  in  Figure  1.  There  are  several  ways  to  advance  the  solution  by  one  A tc  on 
the  grids.  These  are  dependent  on  the  interface  conditions  between  the  grids.  For 
example,  with  a  refinement  ratio  r  =  3,  (r  =  hc/hf  =  the  ratio  of  the  coarse  grid 
Ax  to  the  fine  grid  Ax)  the  integration  order  (from  coarsest  to  finest)  is 
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where  G,  are  grids  at  level  i  with  i  increasing  for  finer  levels.  Grids  at  level  1  are 
integrated  r  times  as  often  as  grids  at  level  i  —  1  but  with  A t,  =  A£j_j/r,  thus  all 
axe  integrated  to  the  same  point  in  time.  The  order  of  integration  assures  that  all 
fine  grids  will  have  sources  for  boundary  values. 

Error  estimation  is  also  repeated  and  solutions  from  the  fine  grids  must  be 
passed  to  the  coarse  grids.  Errors  are  estimated  and  grids  replaced  on  each  level 
after  a  number  of  time  steps  specified  by  the  user.  Grids  at  level  i  will  be  replaced, 
based  on  an  error  estimate  on  level  i  —  1  grids,  r  times  more  often  than  grids  at 
level  i  —  1.  Solutions  on  level  i  axe  updated  with  those  on  a  higher  level  when  the 
higher  level  solutions  have  reached  the  same  point  in  time. 

Errors  on  the  various  grids  axe  estimated  using  a  method  based  on  Richardson 
extrapolation.  If  the  solution  is  smooth  the  local  truncation  error  can  be  expressed 
as 

u(x,  t  +  k)  -  Qh(u(x,t))  =fc(fc?la(x,t)  +  hgib(x,t )) 

+  kO(k9l+1  +  h9i+1) 

=  T+kO(k9'  +  l  +  h1i  +  l), 
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where  91,52  are  the  orders  of  accuracy  in  time  and  space  and  Qh  is  an  operator 
representing  the  finite  difference  scheme  and  defined  as  u(x,t  +  k)  =  Qh(u(x.t)). 
If  we  let  Qih  be  the  same  difference  operator  with  step  sizes  of  2k  and  2 h  and  if 
we  assume  that  the  order  of  accuracy  of  the  time  and  space  differencing  are  equal 
and  also  that  the  solution  is  sufficiently  smooth,  then  we  can  derive  the  following 
expression  for  the  leading  order  term  r  of  the  truncation  error. 

Ql(u(x,t))-QM(u(x,i))  =r  |  0(h<]+2)  (22) 

where  Q\  is  the  operator  Qh  twice  applied  to 

This  gives  us  an  estimate  of  the  local  truncation  error  at  time  t.  The  procedure 
is  to  take  a  giant  step  based  on  mesh  widths  2 h  and  2k  using  the  solution  at  time 
t  and  then  to  compare  it  with  the  solution  found  by  taking  two  regular  steps.  The 
same  solver  used  to  integrate  the  equations  can  be  used  to  estimate  the  error.  The 
estimator  is  independent  of  the  finite  difference  method  and  is  also  indepedent  of 
the  PDE.  The  method  does  not  produce  accurate  estimations  of  r  for  nonsmooth 
solutions  but  no  problem  is  envisaged  because  the  estimates  will  probably  be  large 
in  these  regions  anyway  and  will  lead  to  the  desired  regridding. 

Estimating  the  error  at  grid  points  is  the  first  step  in  the  regridding  procedure. 
The  regridding  procedure  is 

1)  flag  points  needing  refinement, 

2)  cluster  the  flagged  points, 

3)  fit  rectangular  grids  around  the  clustered  points  and 

4)  repeat  if  necessary. 

Gridpoints  are  flagged  if  the  estimated  error  exceeds  a  user  specified  value.  Clus¬ 
tering  the  flagged  points  serves  two  purposes.  First,  it  separates  spaciallv  distinct 
phenomena.  In  many  problems  there  are  often  multiple  shocks  or  fronts.  These 
features  will  then  be  on  different  grids.  The  second  purpose  is  to  subdivide  points 
when  one  large  region  should  be  fit  with  several  grids. 

Clustering  the  grid  points  and  fitting  rectangles  to  the  clusters  are  the  most 
difficult  parts  of  the  regridding  task.  Inexpensive  clustering  algorithms  are  rarely 
satisfactory  for  both  clustering  purposes.  For  this  reason,  a  simple  algorithm  is 
used  for  clustering  in  a  first  pass  and,  if  this  proves  unsatisfactory,  a  more  complex 
and  expensive  algorithm  is  used.  The  simple  method  produces  clusters  according  to 
the  nearest  neighbor  principle.  If  a  point  has  any  other  point  within  some  specified 
minimum  distance  then  these  points  are  in  the  same  cluster.  This  method  works  well 


for  the  first  purpose,  but  very  poorly  for  the  second.  The  more  complex  methods 
use  minimum  spanning  trees  or  nearest  neighbor  graphs.  These  structures  connect 
the  points  in  an  organized  way.  An  iterative  method  is  then  used  which  merger- 
points  with  core  groups  of  points  and  attempts  to  maximize  the  efficiency  of  a  grid 
The  efficiency  of  a  grid  is  a  measure  of  how  large  the  grid  is  compared  to  how  many 
flagged  points  the  grid  contains. 

The  last  task  the  regridding  algorithm  must  complete  is  to  actually  fit  the 
rectangles  to  the  clusters.  There  are  several  methods  which  will  do  this  and  some 
produce,  on  the  average,  more  efficient  rectangles  than  others.  In  the  algorithm  a 
simple,  less  expensive  method  is  used  because  it  works  well  and  also  must  frequently 
be  used  in  the  clustering  routine.  The  method  fits  rectangles  by  computing  a  least  - 
squares  fit  line  to  the  cluster  of  points.  This  line  is  the  principal  axis  of  the  rectangle 
and  an  orthogonal  line  will  be  parallel  to  the  other  axis.  It  is  then  an  easy  matter 
to  compute  where  the  sides  and  the  corners  of  the  rectangles  should  be.  though  this 
is  the  most  expensive  part  of  fitting  the  rectangle.  Finally,  the  rectangle  is  enlarged 
so  as  to  provide  a  buffer  zone  between  the  flagged  points  (the  phenonmena)  and  the 
rectangle  (fine  grid)  boundaries. 

The  data  structure  stores  two  kinds  of  information,  descriptions  of  grids  and 
the  grid  solution  vectors.  It  is  natural  to  think  of  these  grids  in  the  context  of  a 
tree,  with  the  coarse  grid  being  at  the  root  of  the  tree.  At  each  node  of  the  tree  lies 
a  grid.  Each  grid  (node)  is  characterized  by  the  following: 

1)  grid  location, 

2)  grid  point  specifications. 

3)  level  in  tree, 

4)  offspring  pointer, 

5)  sibling  pointer, 

6)  parent  pointer. 

7)  intersection  pointer, 

8)  pointer  to  the  next  grid  at  same  level. 

9)  time  to  which  grid  has  been  integrated. 

10)  index  in  storage  array  for  solution  on  grid. 

This  information  is  stored  for  each  grid  at  the  nodes  of  the  tree.  Figure  2  shows  an 
example  of  the  tree  for  the  grid  system  of  Figure  1. 

All  solution  vectors  are  stored  in  one  array.  This  array  is  managed  as  a  linked 
list  of  used  and  available  blocks  of  storage.  Storage  is  allocated  in  contiguous  blocks 
by  scanning  the  list  of  available  blocks,  taking  the  first  block  that  is  large  enough. 


and  returning  whatever  space  is  unused.  Reclaimed  space  can  be  re-inserted  into  the 
list  and  reused.  The  structure  of  FORTRAN  does  not  allow  for  dynamic  memory 
allocation  outside  the  program,  hence  all  storage  is  defined  initially. 

We  have  described  the  algorithms  which  control  the  integration  sequence,  er¬ 
ror  estimation,  clustering,  gridfitting  and  data  management.  The  program  is  con¬ 
structed  modularly.  The  data  structure  and  the  methods  used  to  alter  it  can  be 
accessed  by  all  routines.  The  user  has  only  to  supply  an  integration  routine  (a 
solver)  which  solves  the  equations  that  he  is  interested  in.  Purposes  necessitating 
changes  usually  require  altering  only  one  or  a  few  modules,  and  not  the  entire  code. 

Our  use  of  the  adaptive  grid  method  and  the  program  just  described  is  greatly 
facilitated  by  code  modularity.  The  program  contains  the  necessary  algorithms  and 
data  structures  to  carry  out  the  adaptive  method  outlined  earlier.  Many  of  these 
algorithms  have  been  borrowed  from  the  fields  of  computer  science,  mathematics 
and  other  disciplines  and  it  is  their  application  for  numerical  weather  prediction 
that  is  new. 
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3.  PRIMITIVE  EQUATIONS 
3.1  Equations 


We  solve  the  Euler  equations  for  dry,  adiabatic  large  scale  atmospheric  flows. 
The  set  commonly  used  for  studying  compressible  fluid  flows  has  as  dependent 
variables  u,v  and  w,  the  horizontal  and  vertical  velocities,  and  T,p  and  p,  the 
temperature,  pressure  and  density,  with  geometric  coordinates  x,  y  and  z.  We 
make  several  changes  to  this  set. 

Large  scale  flows  Eire  very  neatrly  hydrostatic,  and  making  this  approximation 
reduces  the  z  momentum  equation  to  the  hydrostatic  equation  with  the  added 
benefit  of  removing  sound  waves  from  the  solution.  We  also  recast  the  system  using 
the  nondimensional  pressure  o  in  place  the  vertical  coordinate  z.  The  coordinate  o 
is  defined  eis 


a  -p/i: 

where  7r  is  the  surface  pressure.  Thus,  at  the  surface  p  =  ps  =  7r,  o  =  1  and  at 
the  top  of  the  atmosphere  p  =  0  and  a  =  0.  The  vertical  coordinate  o  has  a  range 
0  <  o  <  1. 

Using  the  definition  ir  —  p„  and  a  —  p/ir  we  can  write  the  equation  of  state  as 
on  =  pRT  and  use  it  to  eliminate  the  density  p  from  the  equations.  The  final  step 
is  the  introduction  of  the  variable  <f>  —  gz,  the  geopotential,  which  replaces  z  as  a 
dependent  variable.  The  transformed  Euler  set,  known  as  the  primitive  equations. 


is 


d<f> 

d(lno) 

dn 

dt 


d 


d<t> 


~dx^u  )  “  “  alt™*)  +  *fv  - 


dy 


do' 


dx 


-  RT^-  +  ttFz 
dx 


d(  \  d  (  2n  d,  •  »  f  dd> 

- Tx {l,vu)  -  d-y (TO  >  -  -  *fu  -  Ta? 

-  RT ^  +  ir F, 
dy 


-RT 


Q 

-V „  ■  ( nc„TV )  -  —(nCpTo)  +  n Q  +  RTuj  +  *Ft 
do 

_  do 

—  ■  nV  —  7T- - l-Fir 


(3.1) 


(3.2) 

(3.3) 

(3.4) 

(3.5) 


where 


V  =  ui  +  vj 


11 


R  =  gas  constant 
geopotential  <j>  =  gz 


u>  =  dp/dt  =  tt<7  +  a(dn/dt  +  V  • 

The  independent  variables  sire  x,  y,  a  and  t  and  the  dependent  variables  axe 
and  T.  3.1  and  3.2  are  the  u  and  v  momentum  equations,  3.3  is  the  first  law  of 
thermodynamics  and  3.4  the  hydrostatic  equation.  Equation  3.5  is  the  transformed 
continuity  equation.  A  complete  derivation  of  these  equations  can  be  found  in 
Holton  (1979)  and  Haltiner  and  Williams  (1980). 

The  need  to  make  several  other  assumptions  arises  when  these  equations  are 
used.  The  terms  FX,FV,  and  Q  are  forcing  or  source  terms  that  account  for  processes 
not  explicitly  accounted  for  in  the  dynamics  equations.  The  terms  in  the  momentum 
equations  theoretically  include  the  effects  of  diffusion,  both  turbulent  and  molecular. 
In  the  thermodynamic  equation  Q  represents  latent  and  radiational  heating  and 
cooling,  fluxes  of  heat  from  the  boundaries  and  in  essence  all  diabatic  effects.  Major 
assumptions  underlying  models  often  are  found  in  the  parameterizations  of  these 
terms. 

In  large  scale  flows  the  effects  of  viscosity  and  turbulence  have  negligible  con¬ 
tributions  to  the  forcing  terms  Fx  and  Ft.  Horizontal  diffusion,  fourth  order  in  the 
interior  and  second  order  by  the  grid  boundaries,  is  included  only  to  stabilize  the 
numerical  solution.  This  stabilizing  diffusion  term  is  also  calculated  for  Fv  and  Ft- 
There  is  no  vertical  diffusion  in  any  of  the  equations.  We  sire  solving  for  adiabatic 
flow,  the  model  does  not  consider  radiation  effects  and  there  are  no  sources  of  heat. 

Temperature  inversions  and  instabilities  usually  are  addressed  using  convec¬ 
tive  adjustment  parameterizations.  The  use  of  the  hydrostatic  assumption  in  the 
equations  removes  the  mechanism  which  allows  the  atmosphere  to  respond  to  in¬ 
stabilities  in  a  vertical  air  column.  Convection  (vertical  air  motion)  often  takes 
place  in  the  atmosphere  as  a  response  to  an  instability  in  an  air  column.  The  in¬ 
stability  can  be  thought  of  as  the  presence  of  more  dense  air  above  less  dense  air. 
This  instability  cannot  give  rise  to  vertical  motions  because  there  is  no  feedback  to 
a  vertical  momentum  equation  where  dw/dt  is  driven  by  the  forcing.  A  common 
approach  is  to  parameterize  convection  that  arises  from  these  instabilities  through 
the  use  of  convective  adjustment  schemes.  We  describe  a  simple  dry  convective  ad¬ 
justment  scheme  used  in  this  solver.  For  a  more  detailed  explanation  see  Haltiner 


and  Williams  (1980). 

Air  parcels  can  be  characterized  by  their  lapse  rates  7 


6T 
6z  ' 


If  in  a  dry  atmosphere  the  lapse  rate  at  a  point  is  greater  than  the  dry  adiabatic 
lapse  rate,  then  a  vertical  adiabtic  displacement  of  a  parcel  at  this  point  will  be 
unstable.  This  will  produce  vertical  convection  in  the  region.  A  simple  way  to 
parameterize  this  process  is  as  follows. 

1)  Calculate  the  large  scale  fields  without  considering  instabilities. 

2)  Calculate  the  actual  lapse  rates  and  dry  adiabatic  lapse  rates  in  a  column. 

3)  Compare  the  lapse  rates 

7  <  7 d  stable  6T  =  0 

7  >  7 4  unstable  6T  ^  0 

7<*  is  the  dry  adiabatic  lapse  rate. 

6T  is  temperature  change  in  the  column  due  to  convection 

In  the  first  case  nothing  is  done  for  the  column  is  stable.  In  the  second  case  the 
vertical  temperature  profile  is  adjusted  to  a  neutral  or  slightly  stable  lapse  rate  7 
subject  to  the  condition  that  total  potential  energy  is  conserved,  i.e., 

*Zt 


r£t  c  rPh 

/  cp6Tpdz  =  -M  6Tdp  =  0 
Jzb  9  Jp, 


where  6  and  t  represent  the  bottom  and  top  of  the  unstable  layer  and  cp  is  the 
specific  heat  of  air  at  constant  pressure.  The  scheme  assumes  that  convection  causes 
potential  energy  to  be  converted  to  kinetic  energy  which  is  eventually  converted  into 
internal  energy. 


3.2  Discretization 


On  a  horizontal  plane  the  grid  used  in  this  model  is  the  C  grid  described  in 
Arakawa  and  Lamb  (1977).  The  C  grid  is  shown  in  Figure  3.  It  resolves  shorter 
waves  well,  helps  accurately  represent  wave  group  velocities  and  amplitudes  and 
possesses  very  good  geostrophic  adjustment  properties.  Vertical  differencing  takes 
place  on  the  grid  shown  in  Figure  4.  The  variables  shown  on  the  horizontal  grid 
are  carried  between  a  levels,  a  and  its  time  derivative  are  carried  at  a  levels.  The 
surface  pressure  tt  and  the  surface  geopotential  <f>s  are  defined  at  the  surface.  The 
a  axis  is  always  vertical,  i.e,  radially  outward  from  the  center  of  the  earth,  thus  the 
values  or  gradients  of  tt  or  <p,  at  some  point  on  a  horizontal  plane  above  the  surface 


are  taken  as  those  values  at  the  surface  directly  (vertically)  below  the  point  at  the 
surface. 

We  first  consider  how  the  equations  are  differenced  on  the  horizontal  a  surfaces 
on  the  C  grid.  Centering  a  control  volume  over  a  u  velocity  point  at  (x,y),  where 
x  =  (i  —  1)  •  Ax  and  y  =  (j  —  1)  •  Ay,  we  can  denote  fluxes  of  u-momentum  through 
the  control  volume  faces  in  the  x  and  y  direction  as  Fu  and  Gu .  The  discretization 
for  the  horizontal  advection  terms  in  the  u-momentum  equation  is 


+ £<"»>  -  + 


where 


’  (u«+^.i+^7r«+i.i  +  T<+i,i+i)  +  + 


Vertical  (<r)  derivatives,  for  example  (d(irucr)/d<T),  are  computed  as 
-^-(vucr)  =  - +  u«,>,fc-i)) 

where  the  overline  denotes  an  averaged  value  for  t  and  a  at  the  points  (i,j).  The 
velocity  u  is  averaged  to  arrive  at  a  value  of  u  at  k  ±  | . 

Similarily,  we  can  center  a  control  volume  about  a  n  point  and  denote  “mass” 
fluxes  through  the  surfaces  as  F,  G  and  5.  The  continuity  equation  (3.5)  is  finite 
differenced  as 

d  \  X 

~Ax^'+\ihk  ~  ~  ^i+j-i l,k) 


where 


_  1  ,  . 
*i+±,j,k  ~  2U»+i*>'fc^7r,'>  +  7r*+1’-?' 


Gi,j+},k  ~  2V*J+i,k(7ri,j+l  +7r«,>) 
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The  remaining  terms  in  the  equations  3. 1-3.4  are  differenced  in  a  similar  man¬ 
ner.  The  overall  scheme  will  conserve  mass  (disregarding  the  diffusion  term  in  the 
pressure  tendency  equation)  but  will  not  necessarily  conserve  energy. 

The  leapfrog  method  is  used  to  integrate  the  spacially  differenced  equations 
in  time.  The  method  is  explicit,  second  order  in  time  and  possesses  good  phase  and 
amplitude  wave  propagation  characteristics.  Equations  3. 1,3. 2  and  3.3  axe  marched 
forward  in  time.  The  surface  pressure  ir  is  found  by  integrating  the  continuity  equa¬ 
tion  (3.5)  vertically  at  each  v  gridpoint.  It  is  only  when  integrating  this  equation 
that  the  vertical  boundary  conditions  (<r  =  0  at  <7  =  0, 1)  are  needed.  The  vertical 
integration  of  3.5  results  in 


^  =  £v*-(^)a**. 


The  a's  can  be  found  by  integrating  3.5  down  to  the  required  level.  The  geopotential 
<j>  is  found  by  integrating  the  hydrostatic  equation  3.4  . 

The  stability  of  the  leapfrog  scheme  is  limited  by  the  CFL  condition 


s4*<i 

Ax  ~ 


The  fastest  waves  in  these  equations  are  the  gravity  waves,  where  the  external  grav¬ 
ity  wave  travels  approximately  an  order  of  magnitude  faster  than  the  meteorological 
waves  of  interest.  These  waves  are  important  in  the  geostrophic  adjustment  process 
thus  they  cannot  be  filtered  out  of  the  equations  without  some  adverse  affects.  The 
gravity  waves  severely  limit  the  maximum  time  step  which  can  be  used. 


3.3  Boundary  Conditions 


This  problem  is  posed  as  an  initial  boundary  value  problem  (IVBP).  The  initial 
values  for  the  velocities  u  and  u,  the  surface  pressure  p,  or  tc  and  the  temperature 
T  or  the  geopotential  <f>  are  necessary  as  initial  conditions  for  the  model. 

The  boundary  conditions  which  must  be  specified  are  those  in  the  vertical  (at 
the  top  and  bottom  sigma  layers)  and  in  the  horizontal  at  the  lateral  boundaries. 
The  choice  of  the  sigma  coordinates  leaves  us  with  very  simple  boundary  conditions 
in  the  vertical.  The  conditions  are 


dcr 

cr  =  —  =  0 
dt 


at  both  the  top  ( p  =  0,<7  =  0)  and  the  bottom  (p  =  p3  =  7r,cr  =  1).  At  the  surface  6, 
is  specified  and  this  serves  as  the  lower  boundary  condition  for  the  integration  of  the 


hydrostatic  equation  3.4  .  At  the  upper  and  lower  boundaries  free  slip  conditions 
are  applied  for  the  u  and  v  velocities  and  no  flux  conditions  are  applied  for  the 
temperature. 

The  lateral  boundary  conditions  may  vary  with  the  application  of  the  model. 
Our  test  cases  are  for  flow  in  a  free  slip,  east-west  periodic  channel.  The  north- 
south  boundaries  are  no  flux,  (v  =  0).  We  use  these  boundary  conditions  on  the 
base  (coarse)  grid. 

For  fine  grid  boundary  conditions  we  specify  u,  v,  T  and  n  at  the  boundaries 
using  bilinearly  interpolated  values  from  another  grid.  In  this  regard,  we  are  choos¬ 
ing  to  apply  continuity  conditions  at  the  fine  grid  boundaries  as  opposed  to  treating 
each  fine  grid  as  a  separate  initial  boundary  value  problem. 

3.4  Primitive  Equations  Considerations 

Oliger  and  Sundstrom  (1978)  have  shown  that  the  primitive  equations  are 
ill-posed  for  the  initial  boundary  value  problem  with  open  boundaries.  More  specif¬ 
ically,  they  state  that  “local,  pointwise  boundary  conditions  cannot  yield  a  well- 
posed  problem  for  the  open  boundary  problem  for  the  hydrostatic  equations” .  This 
statement  holds  for  the  continuous  equations  and  is  based  on  examinations  of  the 
behavior  of  appropriate  norms  of  the  solution.  It  logically  follows  that  a  discrete 
approximation  of  these  equations  cannot  have  a  norm  which  behaves  reasonably  if 
it  accurately  approximates  an  ill-posed  problem  (Thomee,  1969). 

Limited  area  modellers  have  traditionally  circumvented  the  problem  of  ill- 
posedness  and  resulting  exponential  error  growth  by  including  horizontal  dissipation 
in  their  models  and,  more  importantly,  by  using  sponge-type  boundary  conditions 
and  increased  horizontal  dissipation  close  to  the  boundaries.  This  leaves  open  the 
question  of  exactly  what  equations  are  being  solved  and  the  accuracy  of  the  limited 
area  model  solutions.  Our  example  calculations  are  for  flow  in  a  periodic  channel. 
The  primitive  equations  are  weakly  well-posed  for  these  boundary  conditions. 

The  question  also  arises  as  to  whether  the  primitive  equations  are  ill-posed  for 
solution  on  the  local  refinements  (fine  grids).  This  would  be  the  case  if  we  chose  to 
treat  a  fine  grid  as  a  seperate  IVBP  and  use  boundary  conditions  appropriate  for 
the  equations.  Instead  we  have  chosen  to  use  continuity  conditions  and  interpolate 
all  data  from  one  grid  on  to  the  boundary  of  another.  Now,  we  must  consider  the 
accuracy  and  stability  of  these  conditions. 

There  tire  no  analyses  for  the  primitive  equations  or  for  nonlinear  hyperbolic 
equations  concerning  the  accuracy  or  stability  of  our  boundary  scheme.  For  a 


2-D  model  hyperbolic  equation  Berger  (1985)  has  shown  that  using  leapfrog  with 
overlapping  grids  and  grid  refinement  in  both  time  and  space  is  stable.  Our  adaptive 
integrations  of  the  primitive  equations  have  also  proven  to  be  stable. 

Accuracy  and  conservation  are  related  issues  and  there  axe  few  results  concern¬ 
ing  overlapping  boundaries  which  axe  rotated  with  respect  to  each  other.  Berger 
(1984)  derives  a  conservative  boundary  scheme  for  use  in  solving  hyperbolic  systems 
of  conservation  laws  on  2-D  rotated  rectangles.  We  do  not  implement  the  scheme 
here  and  know  of  no  implementation  of  it.  Henshaw  (1985)  has  found  that  when 
solving  elliptic  equations  on  overlapping  grids  one  must  use  boundary  value  inter¬ 
polation  schemes  that  are  at  least  as  accurate  as  the  interior  numerical  scheme  and 
in  some  cases  at  least  of  one  order  higher  accuracy.  Browning,  Kreiss  and  Oliger 
( 1973)  show  that  solutions  of  hyperbolic  equations  on  embedded  grids  of  different 
resolutions  may  produce  phenomena  similar  to  that  of  the  propagation  of  waves 
through  materials  of  diferent  densities.  There  can  be  interference  of  waves  which 
have  passed  through  refined  regions  with  those  that  have  not.  This  interference  is 
a  weak  instability  and  obviously  results  in  a  grossly  inaccurate  solution. 

We  address  the  issues  of  accuracy  and  conservation  in  how  and  where  we 
decide  to  place  fine  grids.  In  the  adaptive  scheme,  fine  grid  boundaxy  values  are 
interpolated  bilinearly  from  other  grids.  We  use  an  estimate  of  the  error  in  the 
solution  to  place  the  fine  grids  and  periodically  re-estimate  and  replace  the  fine 
grids  so  that  regions  of  high  error  always  remain  contained  in  the  fine  grids.  The 
regions  of  high  error  (the  disturbances)  must  never  be  allowed  to  pass  through 
fine  grid  boundaries  onto  the  coarse  grid.  Thus  fine  grid  boundaries  always  lie  in 
areas  where  the  solution  error  is  low,  i.e.  where  bilinear  interpolation  will  introduce 
only  small  errors.  Our  integrations  indicate  that  in  this  context  the  use  of  bilinear 
interpolation  is  sufficient  to  insure  accurate  and  stable  solutions. 
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4.  TEST  CASE  RESULTS 
4.1  Barotropic  Cyclone 

The  first  test  case  for  the  adaptive  atmospheric  solver  is  the  simulation  of  a 
barotropic  cyclone  (no  vertical  variation)  which  is  being  advected  by  an  easterly 
flow  in  an  east-west  periodic  channel.  The  equations  are  solved  on  an  /-plane  (/  = 
constant  =  5.0  x  10-4sec-1).  There  is  no  surface  friction  or  energy  sink  other  than 
the  diffusion  used  to  stabilize  the  computations,  hence  the  solution  should  show  the 
cyclone  being  advected  towards  the  west  with  little  change  in  its  appearance. 

The  initial  conditions  axe  shown  in  Figure  5.  Figure  6  is  the  solution  after 
3  days  for  a  coarse  grid  run  with  Ax  =  Ay  =  ISO  km.  and  Figure  7  an  adaptive 
solution  with  a  refinement  ratio  of  three  and  one  level  of  refinement.  These  results 
clearly  show  that  the  coarse  grid  does  not  sufficiently  resolve  the  cyclone  resulting 
in  its  premature  decay  while  the  adaptive  run  does  resolve  the  cyclone  and  shows 
very  little  decay.  The  fine  mesh  is  needed  and  the  adaptive  calculation  is  successful. 

Figure  7  also  shows  the  placement  of  the  fine  grids  in  the  channel.  The  error 
estimation  and  regridding  was  performed  every  nine  hours.  There  is  only  one  fine 
grid  in  the  channel  at  any  one  time  and  the  plots  show  all  the  grids  that  were  placed. 
The  fine  grids  are  usually  just  slightly  larger  than  the  cyclone  and  the  regridding 
occurs  often  enough  so  that  the  fine  grid  tracks  the  cyclone.  Also  note  that  the  fine 
grids  in  this  simulation  are  not  aligned  with  the  coarse  grid.  Another  simulation 
was  performed  where  the  fine  grids  were  aligned  and  had  points  coincident  with 
the  coarse  grid.  The  simulations  differ  only  slightly,  this  demonstrates  that  the 
orientation  of  the  fine  grids  has  little  effect  on  the  solution. 

Fine  grid  placement  depends  on  an  error  estimate  of  the  coarse  grid  solution  as 
described  in  section  2.  Only  the  velocity  error  estimates  were  used  to  place  fine  grids. 
We  will  also  present  error  estimates  for  the  surface  pressure  rr  and  temperature  for 
the  baroclinically  unstable  jet  but  they  are  not  yet  used  for  fine  grid  placement. 
Figure  8  shows  a  typical  error  estimate  for  the  u  and  v  velocity  fields  along  with 
the  fine  grid  placed  over  a  set  of  flagged  points  derived  from  the  error  estimate. 
For  this  flow  the  error  estimate  of  these  fields  places  the  fine  grid  over  the  cyclone 
-  which  is  where  we  expect  it  would  be  necessary. 

Initial  conditions  for  fine  grids  are  interpolated  off  the  coarse  grid  or  fine  grids 
which  existed  before  the  regridding.  A  bicubic  spline  interpolation  is  used  to  obtain 
these  initial  values.  Bilinear  interpolation  has  been  tried,  but  was  found  to  excite 
spurious  gravity  waves  (noise)  which  take  several  hours  to  decay. 


A  critical  component  of  the  adaptive  solution  procedure  is  the  ability  to  accu¬ 
rately  set  the  boundary  conditions  for  the  fine  mesh.  During  this  set  of  simulations 
kinks  arose  in  the  surface  pressure  field  close  to  the  fine  grid  boundaries.  Sur¬ 
face  pressure  errors  are  the  result  of  errors  in  the  mass  divergence  fields  which 
themselves  Eire  the  result  of  errors  in  the  specification  of  the  inflow  velocity.  The 
numerical  scheme  does  not  use  the  pressure  gradient  close  to  the  boundary  and 
the  kinks  do  not  have  any  effect  on  the  solution.  We  include  these  results  only  to 
illustrate  problems  which  can  arise. 

These  errors  in  specifying  the  inflow  velocities  may  have  a  greater  effect  when 
moisture  is  included  in  the  model  and  investigation  reveals  that  the  errors  arise  from 
using  a  numerical  scheme  which  is  inconsistent  close  to  the  boundary.  One  row  in 
from  the  boundary  the  diffusion  is  second  order  as  opposed  to  the  interior  fourth 
order  diffusion.  This  early  version  of  the  code  also  used  a  split  explicit  scheme  for 
advancing  the  gravity  waves,  thereby  allowing  the  use  of  larger  time  steps.  Use  of 
the  scheme  along  with  the  change  in  diffusion  next  to  the  boundaries  was  found  to 
promote  growth  of  the  kink. 

The  adaptive  scheme  attempts  to  minimize  boundary  errors  in  a  manner  not 
connected  to  the  structure  of  the  numerical  scheme.  First  the  fine  grids  are  made 
large  enough  so  that  their  boundaries  are  in  regions  of  small  solution  error,  and 
hence  the  interpolated  values  will  have  small  error.  Second,  if  the  region  of  high 
error  cannot  be  covered  by  a  single  grid  then  multiple,  overlapping  grids  are  used. 
An  important  addition  to  this  is  that  the  fine  grid  boundary  values  must  be  inter¬ 
polated  from  the  best  source,  which  is  often  another  overlapping  fine  grid.  These 
tenets,  along  with  a  consistent  numerical  scheme,  appear  to  be  sufficient  to  produce 
reasonable  solutions. 

Presented  next  are  results  from  the  simulation  of  an  unstable  baroclinic  jet. 
Here  we  use  a  fully  explicit  solver.  By  setting  as  fine  grid  boundary  values  the  first 
two  outer  rows,  we  remove  the  second  order  diffusion  from  the  solution,  and  hence, 
have  a  fully  consistent  numerical  scheme.  Errors  in  the  specification  of  the  inflow 
velocity,  and  hence  the  divergence  and  pressure  fields,  are  greatly  reduced. 

4.2  Baroclinically  Unstable  Jet 

We  have  chosen  for  our  second  test  case  to  simulate  the  evolution  of  an  unstable 
baroclinic  jet  which  has  been  subjected  to  an  initial  perturbation.  The  disturbance 
which  develops  is  commonly  observed  in  the  atmosphere  and  easily  simulated  in  a 
channel.  The  flow’s  close  analog  in  the  atmosphere  and  its  three  dimensional  nature 


allow  testing  of  several  adaptive  code  components  untested  in  the  two  dimensional 
flow  discussed  in  the  previous  section. 

Previous  simulations  of  this  flow  were  performed  in  order  to  gain  an  under¬ 
standing  of  the  basic  physical  processes  driving  it.  Severed  investigators  in  the  late 
sixties  and  early  seventies  (for  example  Williams  1967,  Mudrick  1974)  used  models 
based  on  the  primitive  equations  and  the  quasi-geostrophic  equations  to  simulate  the 
developing  baroclinic  disturbance  and  the  frontal  zones  associated  with  it.  The  cold 
and  warm  fronts  have  received  extensive  analytic  study,  most  notably  in  Hoskins 
and  Bretherton  (1972).  Those  interested  in  the  dynamics  of  this  flow  can  consult 
these  papers  or  for  a  more  recent  and  general  overview  consult  Holton  (19S2). 

For  this  simulation,  we  solve  the  equations  on  a  /3-plane  (Coriolis  parameter 
=  f  =  fo  +  =  df/dy  =constant).  The  grid  has  five  layers  at  a  —  .1,  .3,  .5.  .7 

and  .9.  The  channel  has  a  length  (west-east)  of  5220 km.  and  a  width  (south- north) 
of  S640km.  .  The  north-south  velocity  v  is  initially  zero  and  the  jet  has  no  variation 
in  the  east  west  direction.  The  initially  geostrophically  balanced  jet  is  perturbed 
by  altering  the  north-south  velocity  (u).  After  several  days  simulation  time  a  single 
dominent  wave  appears  in  the  channel.  The  length  of  the  channel  is  then  tripled, 
the  wave  repeated  twice,  and  the  channel  is  now  of  length  L  =  15660km.  .  By 
lengthening  the  channel  we  force  the  adaptive  code  to  use  multiple  overlapping  fine 
grids  -  as  it  might  in  actual  atmospheric  prediction  work. 

The  initial  conditions  for  this  case  are  shown  in  Figures  9-12.  The  jet  core, 
which  contains  the  maximum  jet  velocities,  is  found  on  the  <r  =  0.3  layer.  This 
wave  is  clearly  present  in  in  the  plot  of  the  absolute  vorticity  on  the  a  =  0.5  layer 
(Figure  9).  This  primary  circulation  is  the  result  of  the  baroclinic  instability  which 
arises  from  the  vertical  shear  present  in  the  jet. 

Secondary  circulations  (vertical  and  divergent  motions)  are  driven  by  absolute 
vorticity  advection  and  temperature  advection  in  the  primary  circulation  and  are 
a  result  of  the  hydrostatic  and  geostrophic  nature  of  the  flow.  Positive  vorticity 
advection  occurs  to  the  west  of  the  trough  and  negative  vorticity  advection  to  the 
east.  In  the  lower  layers,  regions  of  cold  and  warm  temperature  advection  are  found 
respectively  at  the  trough  and  crest  of  the  developing  wave.  The  temperature 
advection  can  be  seen  clearly  by  considering  Figure  10.  the  temperature  at  the 
<7  =  0.9  level  and  Figure  11,  the  winds  at  the  same  level.  Horizontal  shear  and 
horizotal  deformation  promote  the  growth  of  the  cold  and  warm  fronts.  The  shear 
and  deformation  fields  intensify  in  the  flow  which  develops  with  the  development  of 
the  surface  pressure  lows  and  highs.  These  surface  pressure  features  are  shown  in 
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Figure  12.  Cyclonic  and  anticyclonic  circulations  form  at  the  lower  levels  around 
the  surface  pressure  lows  and  highes.  These  circulations  are  not  present  in  the  flow 
in  the  upper  layers. 

Poor  representation  of  the  fronts  and/or  of  the  jet  stream  lead  to  slower  devel¬ 
opment  or  even  decay  of  the  disturbance.  The  strength  and  development  of  these 
features  determines  the  adequecy  of  the  grid  resolution.  Figures  13  through  15  show 
the  results  after  three  days  of  simulation  time  starting  from  the  fields  of  Figures 
9-12  and  encompass  three  different  runs;  a  coarse  grid  run  (Ar  =  A y  =  540km.). 
a  fine  grid  run  (Ax  =  A y  =  ISOkm.)  and  an  adaptive  run  with  a  coarse  grid 
Ax  =  A y  =  540 km.,  one  level  of  refinement  and  a  refinement  ratio  of  three  (hence 
A Xfine  «  &Vfine  &  180km.).  Figure  13  shows  the  surface  pressure  for  the  coarse, 
fine  and  adaptive  grid  runs.  On  the  coarse  grid  the  surface  pressure  highes  and  lows 
have  lost  strength  whereas  they  have  not  for  both  the  fine  and  adaptive  runs.  The 
warm  and  cold  fronts  exhibit  similar  behavior.  The  coarse  grid  fronts  are  weak¬ 
ening  while  in  the  fine  and  adaptive  grid  runs  they  are  strengthening.  Again  we 
see  that  the  coarse  grid  cannot  adequately  represent  the  flow  while  both  the  fine 
and  adaptive  calculations  represent  well  the  developing  baroclinic  disturbance.  The 
same  resolution  problem  is  seen  in  the  upper  level  flow.  The  maximum  absolute 
vorticity  associated  with  the  jet  has  grown  from  1.4  x  10_4sec_1  to  1.5  x  10~4sec-1 
after  three  days  for  both  the  fine  and  adaptive  grid  runs  but  it  has  diminished  to 
1.2  x  10_4sec-1  in  the  coarse  grid  run.  The  primary  wave  is  deepening  except  for 
the  coarse  grid  run  where  it  is  being  washed  out. 

Examination  of  the  surface  pressure  fields  in  Figure  13  show  that  the  fine 
grid  run  results  and  the  base  grid  fields  for  the  adaptive  run  results  do  not  match 
exactly.  Indeed  there  are  some  very  noticable  differences,  and  the  differences  are 
even  more  pronounced  in  the  absolute  vorticity  fields.  This  is  simply  because  the 
coarse  mesh  cannot  possibly  represent  all  the  features  that  are  representable  on  the 
fine  mesh.  The  base  grid  fields  for  the  adaptive  run  also  show  the  locations  of  the 
fine  grids  which  have  been  placed  based  on  an  error  estimate  of  the  velocity  fields. 
The  adaptive  fine  grid  fields  match  almost  perfectly  those  of  the  fine  grid  run.  Even 
the  vorticity  fields  (Figure  14)  which  are  sensitive  to  small  errors  in  the  velocity, 
compare  extremely  well  with  the  fine  grid  run  solution. 

In  these  simulations  the  largest  errors  are  associated  with  the  jet.  The  grid¬ 
fitting  routine  fits  just  a  single  grid  over  the  jet  but  this  grid  is  split  into  two 
overlapping  grids  so  that  they  may  be  accomodated  by  the  limited  workspace  in 
the  solver.  The  error  is  re-estimated  every  24  hours  and  new  grids  created  but  the 
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positions  of  the  fine  grids  change  little  over  several  days. 

Here,  as  in  the  previous  simulation,  the  fine  grids  are  not  aligned  with  the 
base  grid.  In  the  overlap  region  the  fine  grids  are  aligned  with  each  other  but  at 
the  periodic  boundary  they  do  not  overlap  in  a  manner  where  their  points  coincide. 
In  both  overlap  regions  solutions  in  the  overlap  region  agree.  For  this  to  be  the 
case  fine  grid  boundary  conditions  must  come  from  the  other  fine  grid  in  regions  of 
overlap.  Using  only  boundary  conditions  interpolated  off  the  coarse  grid  results  in 
unstable  solutions. 

Surface  pressure  contours  run  smoothly  up  to  the  boundary  on  the  fine  grids 
in  the  adaptive  run.  In  our  previous  simulation  there  were  kinks  in  the  surface 
pressure  field  and  large  errors  in  the  divergence  fields  close  to  the  boundary.  By 
using  a  fully  explicit  scheme  (which  allows  the  setting  of  boundary  condition  values 
from  overlapping  grids  in  regions  of  overlap)  and  by  setting  the  variables  at  the  first 
two  interior  rows  as  boundary  values  we  have  eliminated  the  kinks  and  large  errors. 
On  the  fine  grids  the  disturbances  are  well  represented  even  in  the  overlap  regions. 
No  noise  develops  in  the  overlap  region  or  at  the  coarse-fine  grid  boundaries  even 
when  these  boundaries  and  overlaps  have  points  which  are  not  coincident.  Indeed 
it  is  difficult  to  tell  where  the  fine  grid  boundaries  and  regions  of  overlap  are.  Our 
disturbances  remain  as  three  identical  disturbances  even  though  different  parts  of 
the  disturbances  are  represented  in  different  overlap  regions.  Figure  15  is  a  plot 
of  the  temperature  on  the  cr  =  0.9  level.  Here  a  cold  front  and  a  warm  front 
pass  directly  through  fine  grid  overlaps  and  fine  grid  boundaries.  These  fronts  are 
identical  to  those  not  passing  through  an  overlap  region  -  as  they  should  be. 

The  numerical  scheme  used  in  the  solver  does  conserve  mass.  The  adaptive 
method  makes  no  attempt  at  preserving  this  property.  Variables  on  the  base  (coarse) 
mesh  are  updated  where  possible  with  an  averaged  value  from  a  fine  mesh  and  values 
for  boundary  conditions  are  obtained  using  bilinear  interpolation.  For  a  six  day  fine 
grid  run  the  maximum  mass  fluxuation  in  the  channel  is  approximately  .007%  of 
the  total  mass  in  the  channel.  The  adaptive  and  coarse  grid  runs  have  fiucuations 
of  .02%.  The  error  in  the  mass  is  on  the  order  of  a  few  tenths  of  a  percent  at  most, 
i.e.  down  at  the  level  of  the  truncation  error  of  the  numerical  scheme. 

Figure  1G  is  a  plot  of  the  average  kinetic  energy  (I\E)  of  the  atmosphere 
versus  time.  The  integration  is  carried  out  on  the  base  grid  during  the  adaptive 
run  because  the  changing  locations  of  the  fine  grids  make  an  "adaptive"  integration 
difficult.  The  kinetic  energy  should  increase  as  the  disturbances  grow  and  for  the 
fine  and  adaptive  grid  runs  it  does.  Also  note  that  the  oscillations  are  similar  for  the 


fine  and  adaptive  runs.  We  have  calculated  the  KE  which  includes  the  contribution 
from  the  fine  grids  in  the  adaptive  grid  run  at  3  days  and  at  C  days.  These  agree 
well  with  the  fine  grid  values. 

The  total  energy  in  the  channel  (KE+TPE,  TPE  =  gravitational  potential 
energy  and  internal  energy)  is  dominated  by  the  TPE.  Given  that  we  solve  an 
adiabatic  set  of  equations  with  no  energy  sources  or  sinks  and  that  our  flow  is  in  a 
enclosed  channel  we  should  find  that  the  total  energy  is  constant  over  time.  This  is 
not  the  case,  but  in  all  cases  the  gain  in  energy  in  the  systems  is  less  than  a  tenth 
of  a  percent  of  the  average  TPE,  again  on  the  order  of  the  truncation  error  in  the 
scheme. 

One  last  run  was  made  using  two  levels  of  refinement.  Figure  1  shows  the  grid 
arrangment  at  20  hours.  Regridding  at  level  2  was  performed  every  12  hours  and 
level  3  every  4  hours  (refinement  ratio  r  =  3).  The  coarse  grid  has  Ax  =  A ij  — 
540km.,  the  first  level  of  refinement  has  Ax  ss  A y  %  180km.  and  the  second  level 
of  refinement  has  AA'  %  Ay  ss  60km.  Again  the  refinement  is  found  to  be  needed 
around  the  jet  and  the  maximum  errors  are  at  the  jet  core.  A  single  fine  grid 
run  with  Ax  =  60km.,  i.e.  of  the  resolution  of  the  fine  grids  in  the  adaptive  run. 
indicates  that  the  increase  in  resolution  from  ISOkm.  to  60km.  was  unnecessary.  We 
preformed  this  integration  as  a  test  using  multiple  levels  of  refinement.  All  general 
conclusions  found  for  the  two- level  refinement  case  hold  when  using  three  levels. 


5.  ERROR  ESTIMATION 

Fine  grids  are  placed  in  the  solution  domain  based  on  an  estimate  of  the 
truncation  error.  The  procedure  used  to  estimate  the  truncation  error  is  based  on 
Richardson  extrapolation  and  it  has  been  described  in  section  2. 

The  primary  advantage  of  the  Richardson  based  error  estimate  technique  is 
that  the  form  of  the  truncation  error  need  not  be  known.  The  exact  truncation  error 
associated  with  the  discretized  equations  3. 1-3.5  is  complex  and  difficult  to  derive. 
Also  the  leading  order  truncation  error  terms  consists  of  higher  order  derivatives 
which  jure  difficult  to  compute  with  more  than  first  or  second  order  accuracy.  The 
truncation  error  estimate  obtained  using  equation  2.2  with  a  <?th  order  method  is 
accurate  to  0(k1+i  +  h9+1 )  which  for  a  second  order  scheme  such  as  the  one  used  in 
the  present  solver  produces  a  third  order  accurate  estimate  of  the  truncation  error. 

The  error  estimate  for  tne  u  velocity  field  in  the  barotropic  cyclone  case  (Figure 
S)  show  that  the  regions  of  high  error  to  be  around  the  cyclone.  The  coarseness 
of  the  grid  precludes  any  deeper  analysis.  To  further  examine  the  error  estimates 
we  have  estimated  the  truncation  error  for  the  fine  grid  run  at  time  t  =  72  hours. 
The  error  estimates  on  the  fine  grid  contain  detail  which  cannot  be  represented 
on  coarser  grids.  By  comparing  the  Richardson  estimate  with  computations  of  the 
exact  truncation  error  we  can  judge  the  performance  of  the  scheme. 

First  we  should  examine  the  equations  and  directly  estimate  what  the  size  of 
the  truncation  error  should  be.  Our  scheme  is  second  order  in  both  space  and  time. 
Small  time  steps  are  used  in  the  explicit  scheme  due  to  the  presence  of  the  fast 
gravity  waves.  Thus  we  expect  that  the  dominent  truncation  error  will  arise  due 
to  the  spatial  discretization  employed  and  hence  we  will  focus  on  the  error  in  the 
spatial  differencing.  Later  comparison  of  the  spatial  truncation  error  with  the  total 
truncation  error  shows  that  the  spatial  error  does  dominate.  We  can  estimate  the 
size  of  the  truncation  error  by  first  scaling  and  then  nondimensionalizing  equations 
3. 1-3.5  along  with  the  truncation  error.  The  scaling  and  nondimensionalization  of 
equation  3.1.  the  a  momentum  equation,  al  me  with  the  spatial  truncation  error,  is 
contained  in  appendix  1. 

For  large  scale  atmospheric  flows  we  find  that  the  pressure  gradient  and  Cori¬ 
olis  forces  must  balance  each  other  and  that  the  advection  terms  are  an  order  of 
magnitude  smaller  than  these.  This  well  known  result  describes  the  geostrophic 
nature  of  the  atmosphere,  i.e.  the  approximate  balance  between  the  pressure  gra¬ 
dient  and  Coriolis  forces.  Large  scale  flow.-,  can  often  be  considered  in  terms  of 


24 


adjustments  to  maintain  an  approximately  geostrophic  and  hydrostatic  state. 

The  finite  difference  scheme  used  to  discretize  the  equations  is  second  order 
accurate:  the  leading  order  truncation  error  term  is  of  order  Ok(k2  +  h2).  This 
truncation  error  is  the  sum  of  the  truncation  errors  for  the  individual  terms,  all 
of  which  are  of  second  order.  In  the  nondimensional  equations  the  pressure  gra¬ 
dient  and  Coriolis  terms  have  coefficients  of  0(10)  while  the  advection  terms  have 
coefficients  of  0(1).  If  we  look  at  the  order  of  accuracy  of  the  scheme  we  might  ex¬ 
pect  that  the  truncation  errors  for  the  Coriolis  and  pressure  gradient  discretizations 
would  be  an  order  of  magnitude  larger  than  those  of  the  advection  terms.  This  is 
not  the  case. 

The  leading  order  terms  in  the  nondimensional  truncation  errors  for  the  Cori¬ 
olis,  pressure  gradient  and  advection  terms  are: 
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All  primed  variables  are  dimensionless  and  of  0(1)  and  their  derivatives  are  of  0(1); 
thus  all  the  truncation  error  terms  are  of  the  same  relative  size.  Direct  computations 
of  these  terms  confirm  this.  Thus  it  appears  we  cannot  ignore  any  of  the  terms  when 
computing  the  truncation  error  for  the  equations  directly. 

Figure  17  is  a  plot  of  the  truncation  error  in  the  u  velocity  field  at  a  =  .3 
(jet  core)  at  t  =  72  hours  calculated  using  the  Richardson  based  technique.  Figure 
18  is  a  plot  of  the  truncation  error  associated  with  the  spatial  discretization  of 
the  pressure  gradient,  Coriolis  and  advection  terms  computed  using  the  results  in 
appendix  1.  It  is  nondimensionalized  by  multiplying  by  At/(nU),  i.e.,  in  the  same 
way  the  Richardson  estimate  is  nondimensionalized.  The  two  estimates  compare 
very  well.  Both  the  magnitude  and  distribution  of  the  error  are  accurately  predicted 
by  the  Richardson  technique.  We  originally  assumed  that  the  spatial  truncation 
error  dominated  the  overall  truncation  error.  This  comparison  indicates  that  the 
assumption  is  correct. 

We  have  also  estimated  the  error  in  the  surface  pressure  and  temperature 


fields.  The  errors  in  temperature  are  large  where  the  time  rate  of  change  of  the 
temperature  is  large.  This  occurs  in  the  frontal  regions  where  the  temperature 
gradients  and  temperature  advection  are  large.  The  error  estimates  for  the  surface 
pressure  are  very  noisy.  The  noise  may  arise  due  to  imbalanced  initial  conditions 
on  the  2h  error  estimate  grid  which  would  lead  to  gravity  waves  and  hence  noise 
on  the  2h  grid  and  a  noisy  error  estimate.  Also  the  error  in  the  surface  pressure  is 
small  relative  to  the  scaled  surface  pressure  and  the  significance  of  regions  of  high 
error  may  be  small. 


6.  CONCLUSIONS 


An  adaptive  grid  refinement  technique  has  been  used  to  compute  solutions 
to  equations  describing  large  scale  atmospheric  flow.  Fine  grids  are  placed  auto¬ 
matically  based  on  a  Richardson-type  estimate  of  the  local  truncation  error  in  the 
solution.  Simulations  of  a  barotropic  cyclone  and  of  a  baroclinically  unstable  jet 
were  performed  to  demonstrate  the  feasibility  of  using  techniques  of  this  type  in 
NWP  and  similar  large  scale  flow  computations. 

Several  critical  components  ensure  an  accurate,  smooth  solution.  Numerical 
schemes  must  be  consistent  up  to  the  boundaries.  Changing  operators  close  to  the 
boundaries  may  cause  kinks  and  discontinuities.  When  fine  grids  overlap,  boundary 
values  for  one  fine  grid  must  come  from  the  other.  This  necessitates  the  use  of  a  fully 
explicit  scheme  (explicit  even  with  respect  to  the  boundary  conditions)  or  some  new 
scheme  which  would  take  into  account  the  overlapping  fine  mesh  constraint.  Higher 
order  interpolation  techniques  for  use  in  setting  the  initial  conditions  are  necessary 
so  as  not  to  excite  gravity  waves  when  initializing  any  fine  grid.  Given  these  con¬ 
ditions  the  present  solver  conserves  mass  and  energy  for  these  flow  conditions  and 
produces  results  which  compare  well  with  those  on  a  single  fine  grid. 

Richardson  estimates  of  the  truncation  error  compare  well  with  directly  com¬ 
puted  truncation  errors.  The  truncation  error  arising  from  the  spatial  discretization 
dominates  the  truncation  error  and  is  evenly  distributed  among  the  terms  in  the 
equations.  The  error  is  high  in  regions  where  we  would  expect  it  to  be  -  around 
the  cylone  and  the  jet  stream,  and  produces  the  desired  refinement.  No  simpler 
technique  for  estimating  the  truncation  error  is  apparent. 

The  success  of  these  simulations  strongly  supports  the  concept  that  refinement 
should  occur  only  where  dictated  by  the  error  in  the  numerical  solution  and  that 
this  is  sufficient  to  improve  the  accuracy  and  overall  resolution  of  the  entire  solu¬ 
tion.  Using  this  simple  but  key  concept  has  produced  the  first  adaptive  solution  of 
atmospheric  flows  and  the  first  detailed,  quantitative  results  concerning  the  error 
in  the  numerical  solutions. 


APPENDIX  1 


We  wish  to  estimate  the  relative  size  of  the  truncation  errors  associated  with 
the  spatial  discretization.  This  can  be  done  by  scaling  and  nondimensionalizing  the 
equations  and  the  truncation  error  associated  with  the  discretization.  In  this  model 
we  do  not  refine  in  the  vertical  and  our  Richardson  error  estimate  does  not  estimate 
errors  in  the  vertical  differencing.  Thus  to  simplify  the  analysis  we  not  consider  the 
vertical  advection  terms. 

For  our  purposes  it  is  sufficient  to  analyze  just  one  momentum  equation.  We 
scale  and  nondimensionalize  the  u  momentum  equation  (3.1)  with  the  following 
changes  of  variables: 
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where  the  primed  variables  are  dimensionless  and  of  0(1).  The  scaling  values  are 
appropriate  for  large  scale  atmospheric  flows.  If  we  substitute  these  into  equation 
3.1,  drop  the  obviously  lower  order  terms  and  divide  through  by  ~0U/t0  we  arrive 
at  the  following  nondimensional  momentum  equation: 
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Using  these  scalings  we  find  that,  for  for  large  scale  atmospheric  flows,  the  Coriolis 
term  must  balance  the  pressure  gradient  terms  and  the  advection  terms  are  an 
order  of  magnitude  smaller  than  either  of  these.  This  just  describes  the  geostrophic 
nature  of  the  flow. 

We  can  now  use  these  results  to  scale  and  nondimensionalize  the  truncation 
error  associated  with  the  spatial  discretization  of  equation  3.1  .  The  leading  order 


truncation  errors  associated  with  the  terms  in  equation  3.1  are: 
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We  nondimensionalize  the  truncation  errors  by  substituting  the  previously  defined 
primed  variables  and  dividing  by  Un0/t0.  The  leading  order  terms  are  given  in 
section  5.  The  nondimensionalization  of  the  truncation  errors  indicate  that  all  are 
of  the  same  size,  even  though  the  respective  terms  associated  with  the  truncation 


errors  axe  not. 
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Fi|urt  3  Grid  C.  The  variables  T,  4 
and  v  are  found  at  p  points.  «  and  v  are 
the  velocities  in  the  *  and  y  direction 
(east  and  north)  respectively. 


Figure  4  Vertical  structure  of  the  fi¬ 
nite  difference  grid  in  the  sigma  (<r)  co¬ 
ordinate  system  for  a  5  level  model. 
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Figure  5  Initial  conditions  for  the 
barotropic  cyclone.  Surface  pressure 
is  in  millibars.  Winds  travel  counter¬ 
clockwise  around  the  low  pressure.  The 
cyclone  is  advected  towards  the  left  by 
a  zonal  wind. 


Figure  8  Cyclone  after  3  days  with 
the  integration  carried  out  on  the 
coarse  grid  with  Ax  =  180i:m 


Figure  7  Cyclone  after  3  days  using 
adaptive  model.  There  is  only  a  sin¬ 
gle  fine  grid  in  the  region  at  any  time, 
^^ocni  “  180tni.,  — -  60km.  .  Fine 

grids  used  in  the  calculation  are  dis¬ 
played. 


Figure  8  Typical  error  estimate  of 
the  u  velocity  (dimensionless  xlO4). 
The  estimate  has  been  normalized  by 
\U\  mas  =  27  m/ see. 


Figure  9  Absolute  vorticity  (10-Saec-1)  on  the  <r  =  0.5  surface  for  the  baroclinically 
unstable  wave  at  t  =  0.  of  the  adaptive  run. 


Figure  10  Temperature  ( K )  on  the  <r  =  0.9  surface  at  t  =  0.  .  Locations  of  the  warm 
and  cold  fronts  are  shown. 


Figure  11  Velocity  vectors  on  the  a  =  0.9  surface  at  t  =  0.  .  Note  the  positions 
of  the  warm  and  cold  fronts  (shown  in  Figure  10)  and  the  positions  of  the  surface 
pressure  highes  and  lows  (shown  in  Figure  12)  with  respect  to  the  wind  fields. 


Figure  12  Surface  pressure  (millibars)  for  the  baroclinically  unstable  wave  at  t  -  0. 
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face  pressure  in  millibars  at  T  =  72.  hours.  Location  of  the  fine  grids 
:alculation  at  T  =  72  hours  are  shown  in  the  adaptive  coarse  grid 
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Figure  15  Temperature  (/C)  on  the  <r  =  0.9  surface  at  T  =  72.  hours  for  adaptive  fine 
grids  3  (top),  and  the  overlaps  between  fine  grids  2  and  3  (bottom).  The  fine  grids 
overlap  in  the  center  of  the  domain  and  also  at  the  edges  of  the  domain  because  of 
the  periodic  boundary  conditions.  The  solutions  agree  in  the  overlap. 


Figure  16  Horizontally  averaged  kinetic  energy. 


